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ABSTRACT 


A  complete,  unified  description  is  given  ot  the  design, 
implementation  and  use  of  a  family  ot  very  fast  and  etticient 
larye  scale  minimum-cost  (primal  simplex)  network  proyrams. 
The  class  ot  capacitated  generalized  transshipment  problems 
solved  includes  the  capacitated  and  uncapaci tated  generalized 
transportation  problems  and  the  continuous  generalized  assign- 
ment problem,  as  well  as  the  pure  network  flow  models  which  are 
specializations  ot  these  problems.   These  formulations  are  used 
for  a  large  number  of  diverse  applications  to  determine  how  (or 
at  what  rate)  flows  through  the  arcs  of  a  network  can  minimize 
total  shipment  costs.   A  generalized  network  problem  can  also 
be  viewed  as  a  linear  program  with  at  most  two  non-zero  entries 
in  each  column  ot  the  constraint  matrix;  this  property  is  ex- 
ploited in  the  mathematical  presentation  with  special  emphasis 
on  data  structures  for  basis  representation,  basis  manipulation, 
and  pricing  mechanisms.   a  literature  review  accompanies  compu- 
tational testing  ot  promising  ideas,  and  extensive  experimenta- 
tion is  reported  which  has  produced  GfciMNhiT,  an  extremely  effi- 
cient family  of  generalized  network  systems. 


1  .    INTRODUCTION 

This  paper  reports  the  development  ot  a  larye-scaie  primal 
network  code  tor  solving  capacitated  generalized  transshipment 
problems.   The  capacitated  generalized  transshipment  problem  is 
the  most  general  of  the  minimum  cost  tlow  models  in  continuous 
variables,  whicn  include  the  capacitated  and  uncapaci tated  trans- 
portation problems  and  the  continuous  generalized  assignment  problem 
as  well  as  the  pure  network  specializations  ot  these  problems. 
These  models  are  used  tor  a  large  number  ot  diverse  applications 
that  include  transportation  ot  goods,  desiyn  ot  reservoir,  com- 
munications, and  pipeline  systems,  assignment  ot  personnel  ana 
machinery  to  jobs,  bid  evaluation;  currency  exchange  and  cash 
management;  production,  sales  and  inventory  planning;  and  many 
others.   For  further  discussion  of  these  applications  see  survey 
articles  such  as  Bradley  [bj,  or  Glover  et  al.,  [17,  18],  and 
textbooks  (as  well  as  their  cited  references)  such  as  Dantzig 
[13],  Jensen  and  Barnes  [23],  and  Kennington  and  Helgason  [27], 

The  capacitated  generalized  transshipment  model  and  its 

specializations  are  minimum-cost  network  flow  problems.   The 

goal  is  to  determine  how  (or  at  what  rate)  flows  through  the 

arcs  ot  a  network  can  minimize  shipment  costs.   The  network  is  a 

directed  graph,  G  ,  defined  by  a  set  ot  nodes,  N  ,  and  a  set  ot 

arcs,  A  ,  with  ordered  pairs  ot  nodes  (tail,  head)  as  elements 

indexed  by   k:  ( t  , t,  ) .   for  each  arc  there  is  a  shipping  cost 

per  unit  tlow,   c   ,  a  minimum  allowable  tlow  (or  lower  bound), 

I,      ,    and  a  maximum  allowable  tlow  (or  upper  bound,  or  capacity), 

u 
k 


In  addition,  there  are  coefficients  (or  multipliers,  or  gains, 
or  losses),   a    and   b,   which  can  change  the  magnitude  of  (or 

— K  K 

amplify,  or  attenuate)  each  unit  of  flow  respectively  entering 
and  leaving  arc   k  .   Each  node  is  either  a  supply  node  where 
units  of  the  good  enter  the  network,  a  demand  node  where  units 
leave,  or  a  transshipment  node.   The  problem  is  to  minimize 
total  costs  with  flows,   x,  ,  that  satisfy  the  associated  lower 
and  upper  bounds  and  preserve  the  conservation  of  flow  at  each 
node : 


(GN) 


MIN  I      c,  x 


keA 


k  k 


s .  t . 


i 

keA 

with  f,  =i 
k 


a.  x. 
— k  k 


keA 
with  t.  =i 


x,  =  b 


i  e  N 


\  <  xR  <  uR  , 


keA, 


where : 


bi  ■ 


supply  if   i   is  a  supply  node; 
-  demand  it   i   is  a  demand  node; 
U  otherwise. 


Note  that  any  linear  program  with  at  most  two  nonzero  coefficients 
associated  with  each  variable  is  a  generalized  network  (GN).   for- 
mulation (Gin)  is  a  generalized  transshipment  model  (GT)  if  all   a. 
are  +1,  in  which  case  the  corresponding  coefficient  is  called  the 
multiplier   m,  =  -b.   .   For  purposes  of  exposition,  we  will  address 


(GT)  since   assuming  the  existence  of  a  finite  upper  bound  on  each 
variable   it  is  possible  to  transform  the  (GN)  coet t icients--by 
scaling  or  by  reflecting  variables  with  respect  to  their  upper 
bounds--so  that  one  coefficient  is  +1  tor  each  variable.   The  other 
coefficient  for  each  variable  then  becomes  the  multiplier  value, 
m.   ,  and  the  generalized  transshipment  (GT)  network  flow  interpre- 
tation  results  with  a  node  tor  each  constraint  and  a  directed  arc 
tor  each  variable.   If  a  variable  has  two  nonzero  coet t icients ,  its 
arc  is  directed  away  from  the  node  corresponding  to  the  constraint 
with  the  +1  coefficient;  a  variable  with  just  one  nonzero  coefficient 
(±1)  corresponds  to  an  arc  forming  a  self-cycle,  leaving  and 
returning  to  the  same  node. 

Some  generalized  networks  (GN),  such  as  those  obtained  by 
relaxing  integer  restrictions  on  flow  variables,  cannot  be  scaled 
conveniently  to  (GT).   These  (GN)  are  accommodated  by  obvious  minor 
modifications  in  the  following  (GT)  presentation.   We  have  developed 
codes  tor  solving  (GN)  and  specialized  them  tor  ( GT )  problems. 

Generalized  networks  can  be  solved  as  linear  programming 
problems,  but  contemporary  commercial  linear  programming  sys- 
tems consume  much  more  computer  time  and  data  storage  region 
than  special  purpose  network  codes.   Indeed,  the  advantage  of 
network  codes  is  so  pronounced  that  it  is  even  worthwhile  to 
develop  special  linear  programming  procedures  to  exploit  intrin- 
sic network  structure  found  embedded  within  more  general  models 
(e.g.,  Kennington  and  Helgason  [27J,  McBride  [3U],  Brown  and  Graves 
[8],  and  Brown  and  McBride  [9]).   Brown  and  Wright  [11]  and  Brown, 
McBride  and  Wood  [10]  show  that  many  real-life  linear  programs 
contain  a  large  embedded  generalized  network  structure. 
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These  models  are  widely  used  because  they  accurately  de- 
scribe a  large  variety  of  important  applications.   Generalized 
networks  not  only  directly  represent  gains  with   m,  >  1  (e.g. 
interest  return  on  investment,  heat  gain,  etc.)  and  losses  with 
0  <  m,  <  1  (e.g.,  evaporation,  voltage  drop,  attrition,  etc.) 

rK 

in  the  flows,  but  also  admit  conversions  of  units  for  these  flows 
(e.g.,  machine  time  to  output  pieces,  lira  to  yen,  etc.).   In 
addition,   m,  <  0   represents  situations  without  obvious  physical 

K 

flow  interpretation  (e.g.,  flows  which  "enter  the  head"  of  arc   k 
in  proportion   m.   of  those  entering  the  tail),  but  which  nonethe- 
less provide  valuable  modelling  tools.   There  has  been  continuing 
growth  of  interest  in  network  models  because  efficient  computer 
programs  have  made  possible  the  reliable,  economic  solution  of 
problems  with  more  variables  than  virtually  any  other  optimiza- 
tion technique  (e.g.,  the  pure  network  system,  GNET  [6],  has  been 
installed  at  hundreds  of  sites  worldwide  and  is  now  cited  as  a 
routine  research  tool).   Perhaps  most  important,  networks  are 
readily  accepted  by  nonanalysts  and  are  consequently  extremely 
popular  operations  research  models. 

Although  several  papers  have  been  written  in  this  general 
area,  and  significant  computational  breakthroughs  have  been 
reported,  there  has  not  previously  been  a  single,  unified  de- 
scription of  a  complete  implementation,  nor  have  "new  generation" 
computer  programs  been  made  generally  available  to  the  academic 
community.   Here  we  report  the  research  and  computational  experi- 
ments which  have  produced  GfcJNNET,  an  extremely  efficient  family 
of  network  optimization  systems.   GENNET  exploits  pure  network 


structure  embedded  in  generalized  networks,  and  specializes  in- 
trinsically to  GNET  lb],  when  both  systems  are  applied  to  pure 
networks,  (the  floating  point  arithmetic  in)  GENNET  requires 
about  lb  percent  more  time  than  GNET.   An  important  objective  ot 
this  paper  is  to  make  these  new  approaches  easily  accessible  to 
a  wide  audience  via  a  clear  exposition  and  concrete  examples  of 
efficient  FORTRAN  programs.   Further,  the  availability  ot  the 
( GT )  and  (GN)  computer  programs  will  now  make  it  possible  tor 
other  investigators  to  reproduce  and  extend  our  experimental 
results  . 

Bradley,  Brown  and  Graves  lb]  trace  the  historical  develop- 
ments leading  to  contemporary  primal  simplex  pure  network  algo- 
rithms, their  supporting  data  structures,  and  etticient  imple- 
mentations such  as  GNET.   For  other  sub-classes  ot  generalized 
networks,  algorithms  have  been  reported  by  Jewell  [2b],  Eisemann 
[14],  Maurras  [2y],  Glover,  Hultz,  Klingman  and  Stutz  [17,la,ly], 
Balachandran  [2],  and  Jensen  and  Bhaumik  124].   An  efficient 
algorithm  for  large  generalized  network  problems  has  been  de- 
veloped by  Glover,  Klingman,  Hultz,  Stutz,  Karney  and  Elam 
[15,17,18,21,22].   However,  their  contributions  are  scattered 
among  the  papers  referenced;  Kennington  and  Helgason  [27]  and 
Jensen  and  Barnes  [23]  provide  textbook  descriptions  of  compu- 
tations apparently  gleaned  from  these  papers,  providing  an  ade- 
quate treatment  of  the  graphical  algorithm  with  more  computational 
advice  than  the  seminal  presentation  by  Dantzig  [13]  but  tew 
details  ot  etticient  basis  updating. 


2.       THE  APPROACH 

Our  approach  continues  with  generalized  networks  in  pre- 
cisely the  philosophical  vein  of  the  pure  network  exposition 
of  Bradley,  Brown  and  Graves:   we  seek  data  structures  and  algo- 
rithms that  yield  efficient  implementations  without  abandoning 
the  flexibility  of  a  general  large-scale  mathematical  program- 
ming perspective  [6, p.  3  ff].   We  introduce  few  of  the  details  of 
the  general  bounded-variable  simplex  algorithm,  and  we  repeat 
little  of  the  underlying  pure  network  material;  the  assiduous 
reader  might  well  review  the  prior  paper  for  which  this  is  an 
intimate  companion. 

We  continue  with  a  brief  description  of  the  algebraic 
specialization  of  the  simplex  method  for  generalized  networks. 
Specific  design  decisions  and  experiments  carried  out  with  GENNET 
are  described,  including  computational  tests  of  alternate  ap- 
proaches.  Some  extensions  of  GENNET  are  presented  to  further 
exploit  special  problem  structure. 


J.   GENERALIZED  NETWORK  SPECIALIZATION 

Efficient  primal  simplex  specialization  to  the  generalized 
network  case  depends  upon  the  well-known  result  that  any 
generalized  network  basis  can  be  put  in  neirly  (upper)  tri- 
angular form  by  simple  permutation  of  rows  and  columns.   This 
inherent  near  triangularity  can  be  exploited  by  direct  solution 
of  the  simplex  equations  with  modif  ied  forward,  or  back  substi- 
tution.  Fortuitously,  this  basis  structure  also  leads  to  ex- 
tremely fast  solution  updates  orchestrated  in  concert  with 
efficient  dynamic  reorganization  of  each  new  basis. 

Theorem  (e.g.,  Uantzig  [13])  Any  basis  B  extracted  from  a 
generalized  network  problem  can  be  put  in  the  form  (1)  by 
rearranging  rows  and  columns. 


,P 


(1) 


where  each  square  submatrix  component   B   is  either  upper 
triangular  or  nearly  upper  triangular  with  only  one  element 
below  the  diagonal. 


Proof.   Our  proof  is  constructive.   We  know  that  the  basis  has  at 
least  one  nonzero  entry  in  each  row,  and  one,  or  two, 
entries  in  each  column.   Define  a  pairing  as  the  associ- 
ation of  a  row  with  a  column  sharing  an  entry  in   b  ,  a 
deferral  as  a  temporary  deletion  of  a  pair  from  consider- 
ation, and  an  ass  ignment  as  the  fixing  of  a  pair  on  the 
diagonal  in  the  ordering  of  the  rearranged  sequence, 
followed  by  the  assignment  of  all  deterred  pairs  which 
have  a  column  with  an  entry  in  an  assigned  row. 

Step  1)   Defer  singleton  rows.   Locate  a  row  with  one  entry,  pair 
the  row  with  the  column  of  the  entry  and  defer  the  pair. 
Repeat  Step  1  until  no  row  remains  with  one  entry. 

Comment:   Step  1  reduces   B   to  a  submatrix  with  exactly  two 
entries  in  each  row  and  in  each  column.   To  see  this, 
note  that  each  column  can  have  at  most  two  entries,  or 
2m  entries  for   m   columns.   Each  row  has  at  least  two 
entries.   Suppose  that  some  row  has  more  than  two 
entries;  then  at  least  2m  +  1  entries  exist,  leading  to 
a  contradiction.   Thus,  exactly  2m  entries  remain; 
consequently,  each  column  has  exactly  two  entries,  as 
does  each  row. 

Comment:   Step  1  deters  an  upper  triangular  set  of  pairs.   To 

see  this,  diagonalize  the  pairs  in  reverse  of  the  order 

of  their  deferral,  and  place  this  sequence  at  the  end 
of  the  rearrangement  of   B  . 

Comment:   A  sequence  of  assignments  in  the  rearrangement  begins. 


step  2)   Assign  a  Component. 

2.1)  If  a  row  remains  which  is  neither  assigned  nor  deterred, 
begin  a  near  triangular  component:   pair  the  row  and  a 
column  and  ass  ign  this  pair  next  in  the  rearranged 
sequence.   utherwise,  go  to  Step  2.3. 

2.2)  Apply  Step  1  and  ass  ign  each  newly  deterred  pair  next  in 
the  sequence.   Go  to  Step  2.4. 

2.3)  begin  a  triangular  component:   assign  the  most  recently 
deferred  pair  next  in  the  sequence. 

2.4)  Repeat  Step  2  until  all  rows  and  columns  are  assigned. 

Comment:   It  Step  2.1  assigns  a  pair  to  a  component,  then  the 
component  is  nearly  triangular  with  one  element  below 
the  diagonal  in  the  first  column.   utherwise,  the 
component  is  triangular. 

Consider  the  generalized  network  problem  given  in  Figure  1, 

a  generalized  transshipment  problem  with  2  sources,  1U  sinks, 

15  nodes,  and  3U  arcs.   toe  will  use  this  problem  to  illustrate 

concepts  and  etticient  solution  methods.   In  this  problem  the 

multipliers  are  all  positive.   tor  arc   k   directed  trom  node 

i   to  node   j  ,  if  tlow   x,   leaves  node   i  ,  then   ^i^k 

arrives  at  node   j  .   When   m   >  1   the  amount  arriving  at  node 

j   is  greater  than  the  amount  leaving  node   i  .   This  would  be 

the  case  in  cash  flow  problems  when  the  arc  corresponds  to  an 

investment  ana   m,  =  (1  +  r.  )  with   r_,   the  rate  of  return. 

When   m,  <  1   then  the  amount  arriving  at  node   ]   is  less  than 
— k 

the  amount  leaving  node   i  .   here  the  loss  could  correspond  to 


evaporation,  taxation,  transmission  loss,  brokerage  fees,  seepage 

or  deterioration.   Pure  network  arcs  are  indicated  by   m.  =  1 

— k 

(and  may  be  even  more  abundant  in  real  problems  than  in  our 
example,).   for  clarity,  minimum  allowable  flows  are  zero,  and 
maximum  allowable  flows  are  not  specified. 
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1.13 

29 

2 

14 

18.  16 

.98 

3U 

3 

15 

67.62 

1  .00 

Figure  1.   A  Single  Commodity  Generalized  Transshipment 

Problem  ( GT ) 
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figure  2  shows  a  basis  for  the  problem  introduced  in 
Figure  1.   (In  this  simple  case,  there  is  only  one  component: 
p  =  1.) 

A  unique  subgraph  partition  ot   G   denoted   G    corre- 
sponds to   b  .   Let   A   =  {a-^la-1   is  an  arc  associated  with   bJ 

a  column  ot   b}  ,  then   G^  =  [N,A.  J   denotes  the  directed  graph 

£ 
associated  with  the  basis   B  .   To  each  submatrix   b    ot   b 

£       £   £ 
there  corresponds  a  component  ot   G    denoted  by   G   =  [w  ,A  ] . 

£  £        £ 

N    is  the  set  ot  nodes  corresponding  to  the  rows  ot   b    and  A 

£ 
is  the  set  ot  arcs  corresponding  to  the  columns  ot   B   .   It  is 

£ 
known  (e.g.,  [13])  that   G    is  either  a  rooted  tree  or  a  one- 

l 
tree  (a  tree  with  an  additional  arc  forming  one  cycle).   If   G 

£  £ 

is  a  rooted  tree,  then   B   is  upper  triangular.   If   G    is  a 

£ 
is  a  one-tree,  then   B    is  upper  triangular  except  tor  one 

element  below  the  diagonal.   Each  component  can  be  viewed  as 
having  one  cycle  if  we  assert  that  each  rooted  tree  has  a  self- 
cycle  corresponding  to  its  root  node. 

In  our  example  basis,  the  subgraph  G   has  only  one  com- 
ponent ,  (one-tree)  shown  in  Figure  2:   G   is  a  one-tree,  with 
nodes  2,  3  and  11  composing  its  cycle. 

As  in  the  pure  network  case,  this  near  tr iangulat ion  and 
associated  sub-graph  are  naturally  represented  by  a  predecessor 
t_uncticm   p(  ),  and  a  predecessor  graph  (which  does  not  preserve 
the  orientation  of  arcs  in  the  original  network).   The  predeces- 
sor function  can  be  used  iteratively  to  construct  the  unique 
backpath  from  any  node  to  the  root  (or  cycle);  the  backpath 
includes  all  nodes  on  the  cycle.   The  immediate  successors  ot  a 
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Row 
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| 

-1 
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-.79    -1.17 

111                       111 
-1.07 

-1.0 j -.83 

i 

1 

-.91 

-1.0 

-1.0     -1.0 

A  Preorder  Near-Triangulation 
(The  underlined  coefficient  is  below  the  diagonal.) 


Join 


00000       0©0©0 

t  t 

0        0 


One-Tree 


Figure  2.   A  Generalized  Transshipment  Basis 
(For  the  problem  in  Figure  1.) 
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node,  if  any,  are  the  first  nodes  encountered  on  all  paths  ex- 
cept the  backpath  to  the  root,  and  all  the  nodes  on  these  paths 
are  called  successors . 

Note  that  each  basis  may  have  many  near  tr ianyulat ions . 
However,  all  such  near  trianyulat ions  yield  the  same  predecessor 
function  and  graph  (where  the  right  to  left  ordering  of  succes- 
sors of  any  node  is  immaterial).   Thus,  the  predecessor  graph 
does  not  completely  represent  a  near  tr iangulat ion  without 
additional  information:   an  ordering  of  the  rows  (nodes).   For 
algebraic  reasons,  we  restrict  such  partial  ordering  to  preorder 
[6],  in  which  a  node   i   always  precedes  its  successors,  if  any 
and  in  which  all  its  successors,  if  any,  precede  any  node  which 
does  not  precede  node   i  . 
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4.   IMPLEMENTATION 

For  didactic  reasons,  we  bey  in  by  introducing  a  complete 
primal  generalized  network  algorithm  using  a  preorder  traversal 
method.   Controversial  alternatives  are  deferred  until  this 
paradigm  is  presented.   Hereafter,  notation  with  upper  case 
roman  letters  followed  by  parentheses  indicates  a  program  data 
array.   For  instance,  the  predecessor  array  is  referred  to  as 
P(  ). 

Static  arc  storage  is  used  for  tails  T(  ),  heads  H(  ), 
costs  C(  ),  multipliers  MUL(  ),  and  capacities  CP(  ).   Contigu- 
ous storage  by  tail,  or  by  head  node  reduces  T(  ) ,  or  H (  )  to 
an  hierarchical  node-length  entry  point  array.   (GENNET  uses 
contiguous  storage  by  head  node,  as  does  GNET.)   Lower  bounds 
on  arc  flow  are  translated  out  prior  to  solution,  with  appropri- 
ate adjustment  of  the  initial  right-hand  side  of  (GT),  and  of 
CP(  ).   The  sign  bit  of  CP(  )  is  available  to  indicate  arcs 
nonbas  ic  at  their  upper  bounds  ( ref lected  with  flow  -CP(  )  ). 

The  predecessor  function  and  its  array  P(  )  are  defined 
so  that  the  basic  arcs  in  a  cycle  are  oriented  uniformly  in  a 
directed  cycle.   All  basic  arcs  not  on  a  cycle  are  oriented  so 
that  a  backpath  is  created  to  a  cycle.   To  obtain  this  orienta- 
tion the  direction  of  some  arcs  must  be  reversed,  and  the  sign 
bit  of  the  predecessor  array  is  used  to  indicate:   if  P(I)  <  0, 
then  the  orientation  of  arc  (I,  -P(I))  is  reversea  from  its 
original  orientation. 


This  is  the  complement  of  the  discipline  used  in  GNET  [6] 
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A  depth  array  D(  )  reveals  for  each  node  the  number  of 
nodes  on  the  backpath  before  encountering  a  cycle.   Nodes  on  a 
cycle  have  depth  zero.   Number  of  successors,  or  preorder 
distance  are  acceptable  substitutes  tor  depth  [6],  but  are  not 
discussed  here. 

A  preorder  traversal  array  IT(  )  is  maintained  so  that  all 
preorder  successors  of  a  cycle  node  are  encountered  before 
another  cycle  node.   It  is  convenient  to  make  this  a  circular 
list  for  each  near  triangular  basis  component  by  setting  IT(  ) 
of  the  last  preorder  node  in  the  component  equal  to  the  first 
preorder  node  in  the  component. 

The  components  of  G   are  not  inter-connected,  or  equiva- 

i 

lently,  the  sub-matrices  B   in  (1)  do  not  have  common  rows  or 
columns.   Consequently,  the   p   components  of  a  basis  may  be 
represented  in  a  single  set  of  node-length  arrays. 

The  array  X(  )  contains  the  values  of  basic  variables, 
values  of  dual  variables  (or  simplex  multipliers,  or  node  poten- 
tials) are  stored  in  U(  ),  and  IVAR(  )  gives  the  location  of 
basic  variables  in  the  arc  arrays.   The  array  FAC (  )  contains 
the  cycle  factors,  defined  later.   figure  3  shows  these  arrays 
for  the  basis  given  in  figure  2. 

Generalized  networks  do  not  exhibit  totally  unimodular 
bases.   Consequently,  floating  point  representation  is  required 
for  x(  ),  u(  )  and  FAC (  ),  and  is  desirable  for  arc-length  arrays 
C(  )  ,  MUL(  ) ,  and  CP(  ) . 
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Node   Predecessor   Depth   Traversal   basic  Variable   Dual   Cycle  Factor 


P(  ) 

D(  ) 

IT(  ) 

IVAR( ) 

X(  ) 

U(  ) 

FAC  (  ) 

1 

13 

2 

6 

24 

22.860 

16.665 

* 

2 

3 

0 

5 

2 

118.169 

47.148 

-.48101 

3 

11 

0 

1U 

21 

45.826 

31.678 

-.48101 

4 

9 

2 

2 

13 

0.0 

5.418 

* 

5 

-2 

1 

8 

4 

0.0 

27.552 

* 

6 

-3 

1 

15 

7 

21.308 

-3.782 

* 

7 

-2 

1 

11 

10 

3.434 

-38.898 

* 

8 

-2 

1 

14 

12 

27.U87 

27.487 

* 

y 

^3 

1 

4 

14 

9.380 

1.958 

* 

10 

-3 

1 

13 

17 

13.150 

28.606 

* 

n 

-2 

0 

3 

20 

4.170 

-16.053 

-.48101 

12 

-2 

1 

7 

23 

2.076 

32.431 

* 

13 

-3 

1 

1 

26 

11.956 

-13.922 

* 

14 

-2 

1 

12 

29 

22.204 

29.580 

* 

15 

-3 

1 

y 

30 

16.550 

-35.942 

* 

Figure  3.   GENNtT  Basis  Representation  Arrays 
(for  Basis  in  Figure  2) 
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Step  Si,  Priceout 

The  reduced  cost  for  nonbasic  arc  k  ,  oriented  from   f. 

k 

to   t,   is  (given  the  current  dual  solution  u  and  column  N  ): 


r.  =  c,  -  uN 
k     k 


=  ck  "  uf   +  ^kV 
k        k 


=  C(k)  -  U(f  )  +  MUL(k)*U(t,  )  . 

K  K. 


(If  CP(k)  <  0,  arc   k   is  reflected  and  the  siyn  of   r,   is 
reversed.)   At  most,  one  multiplication,  addition  and  subtrac- 
tion are  required.   Note  that  the  multiplication  is  unneccessary 
if  | m,  |  =  1;  further  specialization  is  possible  for  sets  of 
priced  arcs  with  common  attributes.   If  arc   k   is  a  logical  arc 
(slack,  artificial,  or  surplus  variable)  then   C(k)   can  be 
logically  generated,  rather  than  explicitly  stored,  and 


rk  =  C(k)  ±  U(fR 


depending  upon  the  sign  of   P(f  )  . 

K 


i 


From  the  example, 


r2l    =?  C(Z1)    -    U(b)  +  MUL(  zl  )*U(  li) 

=    2U.57  -  (27.552)  +  ,8b(-13.922) 
=  -19.133  . 
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This  variable  will  be  used  as  the  entering  non-basic  variable 
for  further  illustration. 

Step  S2,  Ratio  Test 

For  the  determination  of  the  arc  to  leave  the  basis  the 
system  of  equations 


BZk  =  wk, 


k        k 
must  be  solved  for  the  transformed  column   Z   .   ( -N    is 

used  if  arc   k   is  reflected.)   Due  to  the  near  triangularity 

of   B  ,  this  incoming  column  transformation  can  be  combined 

with  the  ratio  test  in  a  single  integrated  process. 

Suppose  that   N    has  two  nonzero  coefficients  representing 

an  arc  oriented  from   f,   to   t,  ,  and  that  the  arc  is  not 

reflected.   An  apparent  complication  arises  if   f    and   t,   are 

s        t 
in  separate  components  of   B   in  (1),  say   B    and   B   ,  respec- 
tively.  In  this  case  two  disjoint  subsystems  must  be  solved  and 
the  results  added  to  determine  the  nonzero  elements  in   Z 
The  subsystems  are: 


B  (J    =  ef   ,  and 
k 


Bfcg  k  =  -mRe    , 
k 
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fk         tk 

(e    is  the  t,  th  unit  vector;  (J    and  Q    are  disjoint 

k 

k 
components  of  Z   .  ) 

This  complication  is  inconsequential.   In  order  to  see 

this/  consider  solving  one  of  the  systems: 


ut„fck 

b  g   =  -m.e 

k 


The  only  nonzero  elements  of   Q    will  be  those  that  correspond 

to  the  nodes  in  the  backpath  from  node   t   .   As  we  shall  see, 

this  follows  from  the  manner  in  which  the  coefficient  -m,  in 

— k 

row   t    propagates  during  substitution  solution. 

K. 

Suppose  that  G   is  a  rooted  tree.   The  backpath  from  node 
t,   can  be  denoted  by  iteration  of  the  predecessor  function:   t  , 

K.  K 

p(t,  ),  p(p(t  )),...,  root;  this  sequence  is  shown  below  as  a, 

K  K 

d-l,...,0,  analogous  to  the  backpath  length  remaining  to  the  root 


Root 


P(p(tk)) 


ld-2 


d-1 


P(tk) 


ld-l 
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The  values  of   a   and   b   for  each  basic  arc  depend  upon  the 
original  orientation  of  the  arc,  given  by  the  sign  of  P(  ). 
For  original  orientation,  P(  )  >  0  implies  that  a  =  +1  and   b 
is  the  multiplier  value,  P(  )  <  U  reverses  these  definitions. 

The  triangular  system  corresponding  to  this  backpath  is: 


aQ   bx 


ad-2  bd-l 


ad-l  bd 


Q 


U 
0 

• 

0 

u 

-m, 


By  back  substitution,  its  solution  is: 


m, 

— k 

—     d 


g6  =  "b6+lq6+l  = 

36 


Now  suppose  that  G   is  a  one-tree.   Let  node   t,   be  on  the 
cycle.   The  backpath  is   t   ,  p(t.  )  ,  p(p(t,  )),  ...,  c,  with 


p(c) 


t    and  length  p  ;  this  sequence  is  shown  below  as 


-1  ,-2 , . . . , -p ,  analogous  to  the  backpath  length  beyond  the  cycle 
start . 
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The  corresponding  near  triangular  system  is 


a   b   ,, 
-p   -p+1 


'-P  +  1 


a-3  b-2 


~P 


a-2  b-l 


l-l 


0 

0 

0 
0 

-m, 


By  modified  back  substitution,  its  solution  is 


_m-k    1 


"b6+l  q6+l 


for    6  =  -2, . . . ,-p , 


where 


-p   -b. 


f  =  l  -   n 


6=-l  a6 
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The  cycle  factor  (or  loop  factor  [13]),  f  ,  is  common  to  all 
nodes  in  the  cycle  and  can  be  computed  when  the  cycle  is 
created  and  stored  in  FAC (  )  for  all  nodes  on  the  cycle  so 
that  it  is  immediately  available  at  this  step. 

Suppose  that  the  entering  arc  is   ( f  , t.  )   with  coefficient 
entries   (a,,b  )   and  that  the   f,   and   t    backpaths  converye. 
Using  nomenclature  for  the  backpath  sequences  introduced  above, 
the  corresponding  system  is: 


join 


au  bi 


3£-l  b£ 


w 


w 


d-1   d 


ad-X  bd 


0 
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by  back  substitution,  its  solution  is: 


bk 
q,  =  —  ,     ( d  begins  t.   backpath  sequence), 
Q    a  .  k 

—    d 


"b6+lq 


6+1^6+1 


for    6  =  d, d-1 , . . . , r  , 


ld    a 


,     (d  bey  ins  f,   backpath  sequence)  , 


"b6+lq6+ 


tor    6  =  d, d-1 , . . . ,w  , 


<±o- 


i-i 


-(bw^w  +  brV 
a£-l 


b6+lq6+i 
qx  =  tor    6  =  1-2  , 1- i  ,  .  .  .  ,  U  . 


Suppose  that  entering  arc   (f,  ,  t,  )   with  coetticients 
(a,  ,b.  )   enters  and  that  the  backpaths  converge  on  a  cycle.   The 
corresponding  system  is: 
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a-p  b-p+l 
a-p+l 


a-*-l  h-l 


w 


■-* 


a  b 

-s-1     -s 


-s 


a-2  b-l 


~P 


-1 


a     b^i 

w     w+1 

aw+l 


ad-l  bd 


ar   br+l 
ar+l 


0 

\ 
0 


ad-l  bd 
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By  modified  back  substitution,  its  soiution  is: 


q.  =  — —  ,     (d  bey  ins  t,  backpath  sequence)  , 


ld    a 


"b6+i  q6+ 


for    6  =  d, d-i , . . . , r, 


"b  q 


3-8-1 


r^r 

X 


a    ,    t 
-s-1 


"b6  +  i  y6  +  l 


for    6  =  -s-2 ,...,-£, -£-1 ,... ,-s  , 


q-,  =  —  ,  (d  begins  t,  backpath  sequence)  , 
q    a  ,  k 

d 


b6+i  q6+l 


for    6  =  d ,  .  .  .  ,w  , 


<J»-i  * 


-b  q 


*-      a-*-i     t 


D6+l  q&+l  . 

q,.  =  for    6  = 


X-  £  f  %   •   •   f       P  (   ~I-r*«*fO/«*«/   JC   / 


q.  =  q„  +  q,    tor    6  =  -l,...,-p  . 
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For  the  cycle  in  Figure  2, 


i    -(-1)    -(1)    -(-1.17) 
f  =  X J-  x  ^779  x  1 


=  -.481U1  . 

By  now  it  should  be  apparent  that  one  composite  back  sub- 
stitution scheme  will  suffice  for  all  cases.   The  cycle  factor 
is  applied  once  when,  and  if,  a  cycle  is  encountered  on  a 
backpath . 

If  the  backpaths  of   f,   and   t.   converge,  let  join  be  the 
first  node  on  the   t    backpath  that  is  also  on  the   f    backpath. 
If  the  backpaths  converge  on  a  cycle  and  it  the  leaving  arc  pre- 
ceeds  the  join  on  the   f,   backpath,  define  join  as  the  first 
cycle  node  encountered  on  the   f    backpath.   If  the  backpaths  do 
not  converger  join  =  <f>  . 

Several  schemes  are  available  for  identifying  the  join 
efficiently  [6J.   The  depth  (or  number  of  successors,  or  pre- 
order  distance)  of  nodes  on  the  backpaths  can  be  used  to  avoid 
iterating  either  backpath  past  the  join.   Depth,  the  number  of 
nodes  on  the  backpath  until  a  cycle  is  encountered,  can  be  used 
to  indicate  which  backpath  node  is  deeper  and  should  be  iterated. 
When  both  backpath  nodes  have  matching  depths,  zero  depths  indi- 
cate that  each  backpath  is  on  a  root  cycle.   By  remembering  the 
first  root  cycle  node  on  each  backpath,  further  iteration  will 
either  reveal  the  root  cycles  to  be  distinct  with  join  =  $  ,  or 
coincident  with  join  defined  as  above.  When  both  backpath  nodes 
have  matching  depths  greater  than  zero,  the  nodes  are  compared  for 
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equality.   A  match  indicates  the  join,  and  a  mismatch  indicates 
that  both  backpaths  should  be  iterated  tor  another  comparison. 

If  a  join  is  encountered,  the  backpaths  have  merged  and 
either  one  can  be  used  to  complete  the  ratio  test  (GENNET 
continues  the   t    backpath ) .   When  a  cycle  is  encountered  the 
q   values  are  computed  on  the  first  pass  around  the  cycle.   On 
the  second  pass  around  the  cycle  the   q   values  are  computed  and 
the  ratio  test  is  completed. 

As  the  backpaths  are  iterated,  the  column  transformation 
is  applied  and  the  resulting  terms  of  transformed  column,   Z   , 
are  used  in  the  rat  io  test ,  seeking  the  minimum  ratio: 


MIN  i 


CP(k) 


XU) 


'9. 


CP(  IVAKU  )  )-XU  ) 


-  z 


the  capacity  of  the  incoming  arc, 


for    z   >  0,  node  I, 


for    z   <  0  . 


If  a  zero  ratio  is  encountered  during  this  process,  the  ratio 
test  may  be  preemptively  terminated. 

Step  S3,  Pivot 

IF  CP(k)  is  selected  as  the  minimum  ratio,  then  the  enter- 
ing variable  remains  nonbasic  and  is  reflected  to  its  opposite 
bound.   Only  the  flows  X(  )  need  be  updated.   To  do  this,  the 
backpaths  are  iterated  again  ana  for  each  node  I    encountered, 
X(  )  is  reduced  by  z   x  CP(k)  . 
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If  a  basis  exchange  is  required,  an  efficient  update  of 
the  basis  representation  must  preserve  the  rooted  cycle  orien- 
tation, updating  some  entries  in  P(  )  and  D(  ).   Also,  some 
flows  X(  ),  and  some  dual  variables  U(  ),  must  be  changed. 
FAC (  )  must  be  established  for  nodes  on  newly  created  cycles. 
Some  bookkeeping  in  IVAR(  )  and  CP(  )  may  also  be  required. 

The  apparent  intricacy  of  our  task  is  deceiving.   Careful 
analysis  yields  an  elegant  solution.   However,  the  supporting 
arguments  require  close  attention. 

To  simplify  the  explanation,  reorient  the  incoming  arc 
(f,  ,t.)  to  (i,j)  or  (j,i)  (if  necessary)  so  that  the  minimum 
ratio  is  on  the   j   backpath .  Let  the  entering  arc  (i,j)  have 
the  outgoing  arc  (c,d)  on  its  j-backpath.   Also,  reorient  the 
outgoing  arc  (if  necessary)  so  that  the  first  node  encountered 
on  the   j   backpath  is   c  .   Figure  2  shows  a  case  for  which 
both  reorientations  are  necessary. 

In  our  example,  arc  20,  oriented  from  node  2  to  node  11, 
leaves  the  basis  and  arc  27,  oriented  from  node  5  to  node  13 
enters.   We  call  the  backpath  segment  from   j   to  arc  (c,d)  the 
j -stem.   In  Figure  2,  i,  j,  c,  and  d  are  shown.   The  j-stem  is 
composed  of  nodes  5,  2,  3,  and  11.   The  skeletal  update: 

a.  Reverses  the  orientation  of  arcs  within  the 
j-stem;  and 

b.  Urients  the  entering  arc  so  that  it  precedes  this 
redirected  path  with  the  same  orientation. 

If  the   i   and   j   backpaths  merge,  the  node  where  they 
merge  is  called  the  join.   The  join  is  node  3  in  Figure  2. 
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If  the  leaving  arc  lies  beyond  the  join  on  the  backpaths 
a  new  cycle  is  created  in  the  basis  exchange.   In  this  case 
the  portion  of  the   i   backpath  from  node   i   to  the  node  with 
the  join  as  its  predecessor  is  called  the  i-stem.   When  node   i 
is  on  the   j   backpath  the  i-stem  is  null.   In  Figure  2,  the 
i-stem  is  node  13. 

Figure  4  displays  the  pivot  logic  to  be  applied.   The 
algorithm  visits  each  node  affected  by  the  basis  update 
exactly  once .   It  proceeds  up  the  j-stem  one  node  at  a  time 
visiting  the  preorder  successors  of  each  stem  node  via  IT(  ). 
If  a  join  is  encountered  it  switches  to  proceed  up  the  i-stem 
one  node  at  a  time  and  then  returns  to  the  j-stem.   At  each 
j-stem  node,  the  successors  of  the  next  lower  stem  node  have 
already  been  visited.   The  unvisited  successors  of  the  current 
stem  node  can  be  divided  into  two  groups:   the  left  successors 
are  the  nodes  visited  in  preorder  by  iterating  IT(  )  from  the 
current  stem  node  until  the  next  lower  stem  node  is  encountered, 
and  the  right  successors  of  the  stem  node  are  the  remaining 
unvisited  nodes  reached  by  further  iteration  of  IT(  ).   In  Figure 
2,  nodes  8,  14,  12,  and  7  are  right  successors  of  node  2  and 
nodes  10,  13,  1,  6,  15,  9,  and  4  are  left  successors  of  node  3. 
As  we  climb  the  i-stem  the  traversal  IT(  )  is  moditied  so  that 
the  last  of  the  left  successors  (if  any)  points  to  the  first  of 
the  right  successors  and  the  last  of  the  right  successors  (it 
any)  points  to  the  previous  stem  node.   As  we  climb  the  j-stem, 
the  traversal  is  modified  so  that  the  last  of  the  left  succes- 
sors (if  any)  points  to  the  first  of  the  right  successors  and 
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f       BITEH         J 


Set  first 

J-stem  node 

to  visit. 


Visit  left  successor 


Figure  4.     Pivot  Traversal   Scheme 
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the  last  of  the  right  successors  (if  any)  points  to  the  next 
node  up  the  stem  (because  the  update  reverses  predecessor 
orientation  for  the  j-stem). 

The  immediate  switch  to  the  i-stem  upon  encountering  the 
join  during  the  j-stem  iteration  is  motivated  by  a  subtle 
complication:   while  visiting  the  left  and  right  successors 
of  the  j-stem  nodes,  the  nodes  on  the  i-stem  and  their  succes- 
ors  must  be  skipped  if  encountered.   Because  the  i-stem  nodes 
are  successors  of  the  join,  visiting  the  i-stem  as  soon  as  the 
join  is  encountered  (if  one  exists)  on  the  j-stem  leaves  us 
with  the  preorder  successor  of  the  last  i-stem  node  visited. 
This  valuable  artifact  enables  the  subsequent  j-stem  iteration 
to  immediately  skip  all  i-stem  nodes  and  their  successors  should 
they  be  encountered.   This  is  the  key  step  preserving  an  effi- 
cient one-pass  basis  update.   In  Figure  2,  i-stem  node  13  and 
its  successor  node  1  are  successors  of  the  join,  node  3.   The 
preorder  successor  of  the  last  i-stem  node  (called  the  preorder 
link  in  Figure  4)  is  node  b. 

A  stem  node  may  have  a  right  successor  which  is  on  the 
root  cycle  (with  depth  zero).   The  preorder  traversal  array  is 
organized  so  that  all  successors  of  a  cycle  node  are  encountered 
before  the  next  cycle  node.   This  implies  that  a  cycle  being 
broken  by  a  leaving  arc  will  always  be  encountered  as  a  right 
successor  of  a  stem  node.   In  Figure  2,  the  broken  cycle  is 
encountered  as  a  right  successor  of  node  2;  if  arc  (2,11)  were 
not  the  leaving  arc,  then  node  11  would  be  a  preorder  successor 
of  node  2. 
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The  basic  arc  flows,  X(  ),  are  changed  as  each  arc  is 
visited  on  the  j-stem,  and  (if  a  new  cycle  is  created)  on  the 
i-stem.   If  no  cycle  is  created,  X(  )  is  changed  only  it  the 
minimum  ratio  is  nonzero,  and  then  only  on  the  arcs  visited  on 
the  backpaths. 

During  the  update,  the  dual  variables  must  be  changed  so 
that 


ck  =  u    -  mR  u    , 
k         k 


for  every  basic  arc   k   oriented  from  node   t    to  node   f   , 

K  K 

With  incoming  arc  (i,j)  (reoriented  as  in  Figure  2  so  that  the 
outgoing  arc  is  on  the  j-stem),  this  relationship  is  retained 
for  all  nodes  except  for  those  which  the  update  changes  to  be 
successors  of    i:   (i.e.  the  nodes  on  the  j-stem  and  their 
successors ) . 

If  a  cycle  is  not  created,  this  update  proceeds  for  each 
j-stem  node  and  its  successors  as  these  nodes  are  visited  in 
preorder.   for  node   s  ,  and  associated  basic  arc  £  oriented 
from   s   to  p (s )  , 

U(s)  =  CU)  +  MULU)  *  U(P(s)  )  , 

while  the  reverse  orientation  -P(s)  to   s   yields 

U(s)  =  (CU)  -  U(-P(s)))/(-MULU))  . 
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As  in  pure  networks,  the  preorder  traversal  assures  that  a  value 
of  the  dual  variable  of  a  predecessor  node  is  always  determined 
prior  to  its  use  by  any  immediate  successor  nodes. 

When  a  cycle  is  created,  the  dual  variable  must  be  deter- 
mined for  one  of  the  cycle  nodes.   Then,  the  dual  variables  of 
the  remaining  cycle  nodes  and  their  descendents  can  be  found 
one  at  a  time  in  preorder  traversal. 

This  key  dual  variable  is  computed  immediately  after  the 
ratio  test  predicts  creation  of  a  new  cycle  (the  leaving  arc 
lies  beyond  the  join  on  the  ratio  backpaths).   Consider  the 
current  context  for  a  new  cycle: 


rr 


/        * 


P(i) 


P(j) 


from  which  the  new  cycle  will  be  formed  (with  modified  pre- 
decessor function  p  (  ). )  : 


/a  p<p(j>> 

rv rx 


p(j) 


^ — ^d     , .,  a  „v        b  ^  a  ,,  ^        b  , 


-p+1  d-3  b-2  a-2 


-  P 


-1 
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The  corresponding  near-trianyular  system  is: 


(u-p'u-P+i 


/u_1) 


a    b    , 
-p   -p  +  1 


-p  +  1 


-p 


a_3  b_2 


a-2  b-l 


=  (c   ,c 

-p   -p  + 


i.1r.../C_i) 


where  the  indexing  of   c   is  understood  to  yield  basis  arc 
costs  (accomplished  by  indexing  with  IVAR(  )  ). 

The  determination  of  one  term  (say,  u_, )  of  the  solution 
of  this  system  is  induced  by  modified  forward  substitution  to 
be  (e .g  .  ,     [ 11]  )  : 


-      U'  ,    X   —   , 


-1       ""I 


where  (the  cycle  factor)   f   is  defined  (again)  as 


-b 


f  =  1 


6=-l 


and  u'   is  the  corresponding  term  of  the  solution  of  the  strictly 

upper  triangular  system  (omitting  the  cyclic  coefficient  b    , 

-p 

and  using  forward  substitution): 
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c 

-p    a-p 


i     &  6   6-1  r      ,  ,      , 

u;  = for     6  =  -p+l,...,l  . 

o         a  „ 


Note  that  computation  of  u'   requires  that  we  traverse  the 
new  cycle  in  a  direction  opposite  to  its  predecessor  orientation. 
However,  before  the  update  creating  the  new  cycle,  the  j-stem 
exhibits  proper  orientation  for  at  least  part  of  our  work.   Thus, 
we  can  complete  the  first  portion  of  the  forward  substitution 
for  u\    and  accumulate  the  associated  partial  product  component 
of   f   while  iterating  the  j-stem  before  the  update. 

The  remainder  of  the  new  cycle  is  accessed  by  iterating 
the  i-stem.   As  we  proceed  up  the  i-stem  the  remaining  product 
terms  of   f   are  accumulated,  and  the  reverse  i-^tem  path  is 
stored  (e.g.,  using  U(  )  locations,  which  contain  obsolete  dual 
values  to  be  replaced  during  the  imminent  update).   Reaching 
the  join,  this  stored  reverse  i-stem  path  is  then  accessed  to 
complete  computation  of  u'   . 

The  newly  created  cycle  in  Figure  5  is  composed  of  j-stem 
nodes  5,  2,  and  3,  and  i-stem  node  13.   The  dual  system  becomes 


~u13  +  U3  =  4b,6(J 


u   +  u  =  lb  .47 


u   -  -74u5  =  26.76 


88u  +     u   =  20.67  , 

13  5 
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which  yields 


uj3  =  -19.0142, 


f  =  0,3488  (the  new  cycle  factor),  and 


u,~  -  -54.5132  (the  new  dual  for  node  13) 


Thus,  when  a  new  cycle  is  to  be  formed,  the  new  cycle  nodes 
must  be  visited  once  (after  the  ratio  test  and  prior  to  the 
pivotal  update). 

The  one-pass  preorder  traversal  update  can  now  proceed  as 
presented  in  Figure  4.   The  basis  representation  arrays  are  all 
modified  on-the-fly  during  this  traversal.   The  update  of  nodes 
on  a  newly  created  cycle  (if  any)  includes  establishing  the  new 
cycle  factor  PAC (  ).   Changes  for  P(  ),  D(  ),  IT(  ),  IVAR(  ), 
X(  ),  and  U (  )  proceed  analogously  to  the  pure  network  case 
(e.g.,  GNET  [6])  with  simple  modification  of  generators  for  X(  ) 
and  U(  )  to  accommodate  generalized  network  coefficients  for 
basic  arcs. 

Figure  5  shows  the  new  basis  (derived  from  the  example  in 
Figure  2)  before  restoring  near-triangulat ion  with  the  update. 
A  new  cycle  is  to  be  created  and  the  new  cycle  factor  and  a  dual 
solution  (for  node  13  on  the  i-stem)  are  found  at  this  stage. 

Figure  6  displays  the  new  basis  restored  to  near-triangular 
form.   At  this  point,  all  basis  representation  arrays  are  updated 
with  the  values  shown  in  Figure  7  (data  in  italics  has  been 
changed  by  the  update  )  . 
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Figure  5. 


New  Basis  Before  Restoring  Near-Triangulation 
(entering  column  27,  arc  (5,13)) 
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Figure  6.   New  Basis  Near-Triangulated 
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GENNET  is  designed  to  exploit  intrinsic  pure  network  struc- 
ture commonly  found  embedded  in  generalized  network  problems. 
Note  that  when  this  algorithm  is  applied  to  pure  network  prob- 
lems, _i_t  automatically  adapts  to  a  minor  variant  ot  GNET.   of 
course/  GENNET  uses  floating  point  arithmetic  operations  which 
are  intrinsically  slower  on  most  computers  than  the  pure  additive 
integer  arithmetic  of  GNET  (also,  floating  point  arithmetic 
requires  some  extra  editing  for  mantissa  truncation  errors). 

To  mitigate  this  disadvantage,  GENNET  can  test  logically  for 
pure  network  arcs  (with  unit  multipliers)  and  avoid  unnecessary 
floating  point  multiplication  and  division  operations. 

GENNET  also  employs  an  automatic  dual  basis  aggregation 
refinement  ([6], p. 26  f f ) .   Explicit  values  for  D(  ),  U(  )  and 
IT(  )  are  maintained  only  for  nodes  with  successors.   An  array, 
A(  ),  records  for  each  node  the  current  number  of  its  aggregated 
successors .   When  an  aggregated  node  is  encountered  in  the 
priceout,  its  dual  is  generated  from  that  of  the  immediate 
predecessor  of  the  node.   When  a  backpath  ot  an  entering  arc 
begins  with  an  aggregated  node,  it  is  disaggregated,  and  when  the 
leaving  arc  isolates  nodes  with  no  successors,  they  are 
aggregated . 

Figure  8  shows  the  arrays  affected  by  the  dual  aggregation 
scheme  for  the  basis  in  Figures  2    and  3.   IT(  )  indicates  aggre- 
gated nodes  with  (J  entries,  and  these  nodes  have  broken  outlines 
in  the  one-tree  depiction.   The  priceout  of  arc  (5,1J)  requires 
that  U(5)  be  generated  using  the  predecessor  dual  U(2).   Node  b 
subsequently  starts  the  j-backpath  and  is  disaggregated.   The 
update  leaves  node  11  with  no  successors,  and  thus  aggregated. 
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Figure  7.   GENNET  Basis  Representation  Arrays 
(for  Basis  in  Figure  5) 
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Figure  8.   Aggregated  Basis  Representation 
(for  Figures  2  and  3) 
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b.   COMPUTATIONAL  EXPERIENCE 

Significant  design  alternatives  tor  GENNET  have  each  been 
evaluated  by  extensive  experimentation  at  large  scale.   Illus- 
trative computational  experience  is  abstracted  in  this  section 
for  some  of  the  prototype  systems  tested.   Departing  somewhat 
from  the  style  of  the  paper  documenting  such  work  for  GNET  lb], 
relative  performance  is  reported  even  lor  some  competitive  design 
features  subsequently  rejected  tor  adoption  (some  readers  of  [bj 
have  concluded,  quite  incorrectly,  that  only  those  teatures 
reported  tor  GNET  were  tested).   This  should  help  other  re- 
searchers avoid  our  mistakes,  and  may  even  change  some  widely 
held  misconceptions  and  correct  a  few  translation  errors  in 
textbooks  . 

Among  the  key  issues  to  be  resolved  are: 

a)  Static  Storage  of  Arcs .   The  arc  lists  can  be  stored 
in  arbitrary  sequence,  or,  to  save  space,  arcs  can 
be  stored  contiguously  by  tail  node,  or  by  head 
node,  thereby  replacing  an  arc-length  index  array  by 
a  (presumably  much  shorter)  hierarchical  node-length 
entry  point  array. 

b)  Preorder  Manipulation  of  Basis .   The  triple  label  (aug- 
mented predecessor  index)  method  [zu,22\    will  be 
presented  and  compared  with  the  preorder  traversal 
method . 

c)  basis  Aggregation.   An  aggregated  basis  representation 
will  compete  with  an  explicit  representation. 
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d )  i^lciny  Schemes .   Candidate  list  schemes  and  explicit 
arc  pricing  mechanisms  widely  used  in  general  linear 
programming  systems  will  vie  with  dynamic  candidate 
queue  disciplines. 

e)  Pure  Network  Specialization .   Generalized  network 
algorithms  would  ideally  adapt  to  pure  networks  with 
efficiency  comparable  to  pure  network  codes. 

t ■)   Starting ,  Tuning  and  Tai lor ing .   Which  algorithm 

parameters  and  settings  lead  to  high  efficiency  tor 
interesting  classes  of  problems?   Are  heuristics  tor 
advanced  starting  solutions  worthwhile? 
^ )   Generalizations.   Advanced  teatures  and  generalizations 

will  be  suggested. 
Computational  tests  have  been  made  with  many  problems, 
including  a  benchmark  suite  ot  pure  network  problems  generated 
with  NETGEN  [28J,  and  generalized  network  problems  generated 
by  NETGENG  [17,1»J.   figure  y  gives  some  problem  characteristics 
Static  arc  storage  has  been  implemented  in  three  ways: 
Contiguous  by  head  node  with  hierarchical  node- 
length  entry  point  array, 
Contiguous  by  tail  node,  with  hierarchical  node- 
length  entry  point  array,  and 
Explicit  arcs  in  arbitrary  order. 
In  competition  with  our  basis  manipulation  using  preorder 
traversal,  a  triple  label  representation  originally  suggested 
by  h'.  Johnson  [2bJ  has  been  implemented  tor  pure  networks  and 
called  the  augmented  predecessor  indexing  method  by  Glover, 
Karney,  and  Klingman  [2UJ.   Glover,  Klingman  ana  stutz  122J 
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NETGEN 
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NG15 
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200 

200 
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20 
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8 
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20 
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40 
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8 

60 
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40 
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4 
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80 
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4 

12 
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80 
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50 

50 
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0 
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50 
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0 
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50 
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bO 
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0 
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75 

75 
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0 
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75 

75 
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1,500 

75 

75 
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0 
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75 

75 
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0 
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200 
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100 
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0 
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100 

100 
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100 

GTU7 
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115 

4,000 

0 

GT12 

400 

20 

100 

5,000 

0 

GT  1  5 

1,000 

5 

yy5 

4,000 

100 

GTlb 

1  ,000 

20 

100 

6,000 

100 

GT  1 8 

1,000 

30 

400 

7,000 

0 

Figure  y.   bpme  Benchmark  Problems 

(Problem  NG27  is  extensively  studied  in  [6J 
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report  that  the  method  has  been  extended  to  generalized  networks, 
but  reveal  no  details. 

We  have  implemented  our  own  elticient  version  of  the  triple 
label  scheme. 

The  triple  label  representation  uses  predecessor,  successor, 
and  brother  pointers  lor  each  node.  Figure  1U  shows  these  arrays 
lor  the  basis  in  Figure  2. 

To  briefly  illustrate  the  triple  label  scheme  tor  gener- 
alized networks,  let  our  situation  be  the  same  as  tor  the  pre- 
order  example  (the  j-backpath  of  the  incoming  arc  is  arranged 
to  include  the  outgoing  arc  and  to  encounter  node   c   first  on 
that  backpath).   Let 


Vd+1  =  X  ' 


and  the  backpath  from   j   to   c 


vd'  Vd-1 C 


The  skeletal  update  scheme  is 
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6 

-3 

o 

lb 

7 

-2 

u 

0 

8 

-2 

u 

14 

9 

-3 

4 

0 

10 

-3 

U 

13 

11 

-2 

3 

b 

12 

-2 

U 

7 

13 

-3 

1 

b 

14 

-2 

u 

12 

15 

-3 

u 

9 

Figure  lu.   Triple  Label  Basis  Representation 
(  for  t'  igure  2  ) 
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1.  Set  6  =  d 

2.  It  s(v6_^)  =  v6'  ^oto  Step  J 

Otherwise,  it  possible,  tind  a  node  v* 

such  that  P(v*)  =  v 


6-1 


and   B( v* )  =  vr , 

— — > —  o 


and  set  B(v*)  «■  b(v  ),  Go  to  Step  4 

3.  Set  b(v(5_1)  +  B(v6  )  , 

4.  Set  P(v6)  ♦  v6+1, 

B<V  *  S(v6+1}' 

b.   Set  6  <■  6-1, 

Then,  it  vx  +  i  *  c,  Go  to  ^tep_  2, 
Otherwise,  Stop. 

(Step  2  exhibits  the  key  extension  tor  generalized  networks  of 
the  pure  network  triple  label  scheme.) 

These  tive  steps  reverse  the  orientation  ot  all  arcs  on 
the  j-stem  and  orient  (v^+_i'vh)  so  that  it  beyins  this  redirected 
path.   All  other  triple  label  operations  are  obvious  alterations 
of  the  preorder  traversal  procedure. 

A  static  candidate  list  pricing  strateyy  (e.g.,  [I7,la,jl])) 
an  explicit  arc  pricing  method  reminiscent  ot  yeneral  linear 
proyraioming  systems,  and  a  dynamic  candidate  queue  lb]  have  been 
tested. 

The  (L,,L.J  candidate  list  procedure  is  a  simple  strategy. 
Nodes  with  leaving  arcs  (or  entering  arcs  for  contiguous  head 
node  arc  lists)  are  sequentially  priced,  placing  the  most  nega- 
tive candidate  (it  any)  trom  each  node  on  the  candidate  list 
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until  L,.  candidates  have  been  located  (arbitrarily  organized  arc 
lists  require  explicit  arc  pricing).   Entering  arcs  are  chosen 
trom  the  candidate  list  by  most  negative  reduced  cost  until   L, 
iterations  have  been  carried  out/  or  until  all  list  entries  have 
non-negative  reduced  costs.   Arcs  with  non-negative  reduced  costs 
are  dropped  from  the  list  and  replaced  by  continuing  to  scan  the 
nodes  until  L..    candidates  have  been  tound.   It  an  exhaustive  pass 
through  the  nodes  results  in  less  than  L   candidates,  then  an 
optional  closing  gambit  sets  L   equal  to  the  actual  candidates 
found/  and  reduces  L,  by  halt  unless  L,  equals  one.   Glover/ 
Hultz,  Klingman  and  Stutz  [ 1 7 , 1 8  J  report  (L,,Ly)  ot  (b,lU)  to  be 
best  in  their  work. 

The  candidate  queue  is  a  dynamic  list  ot  interesting  arcs 
and  nodes,  scanned  in  a  cyclic  manner.   The  entering  arc  is 
selected  from  the  queue  by  pricing  nne  entries;  if  an  interesting 
node  is  encountered  it  is  replaced  by  its  best-priced  entering 
arc  (or  leaving  arc  tor  contiguous  tail  node  arc  lists).   Arcs 
pricing  favorably  are  retained  in  the  queue.   When  the  end  of 
queue  is  encountered/  the  queue  is  refreshed  by  pricing  1PG  nodes 
in  a  cyclic  general  arc  scan.   During  an  opening  gambit  ot  nns 
pivots,  the  nodes  incident  to  the  entering  basic  arc  are  added  to 
the  queue.   There  is  no  closing  gambit,  since  the  queue  automa- 
tically shrinks  and  finally  collapses  at  optimality.   Bradley, 
brown,  and  Graves  16J  suggest  WNhi  =  62,    NNh  =  3m/4  and 
1PG  =  m/lU  +  1  for  pure  network  problems  with   m   nodes. 

A  rule  to  break  ties  in  the  ratio  test  which  guarantees 
finite  convergence  for  pure  transshipment  problems  has  been 
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developed  by  Cunningham  [12].   Bradley,  Brown  and  Graves  16] 
show  that  the  conditions  necessary  for  finite  convergence  are 
naturally  satisfied  by  GNET  on  over  yu%  of  its  degenerate  pivots. 
Elam,  Glover  and  Klingman  lib]  have  observed  that  the  results  of 
Cunningham  can  only  be  extended  to  the  generalized  network  case 
when  the  multipliers  are  positive.   We  have  not  used  Cunningham's 
modif  ication . 

Bland  [4]  presents  a  class  ot  restrictions  of  pricing  and 
ratio  tests  for  general  linear  programs  which  relies  exclusively 
on  primal  simplex  representation  and  guarantees  finite  conver- 
gence.  These  rules  are  easily  modified  to  produce  an  efticient 
finite  simplex  algorithm.   The  modification  interferes  with  ef- 
fective pricing  strategies  only  during  degenerate  pivot  sequences, 
and  the  restrictions  increase  in  severity  only  with  the  number  of 
pivots  in  that  sequence.   However,  during  a  degenerate  pivot  se- 
quence restriction  records  must  be  accumulated  (e.g.,  a  list  with 
each  incoming  variable  in  one  ot  our  schemes ) .   This  record  is 
naturally  accommodated  by  the  dynamic  candidate  queue,  but  not  by 
a  static  candidate  list  or  explicit  pricing.  No   purpose  is 
served  by  reporting  such  unbalanced  competition. 

A  starting  strategy  has  also  been  tested  in  conjunction  with 
pricing  alternatives.   A  straightforward  starting  method  ex- 
amines each  node  with  supply  (or  demand  tor  contiguous  arc  stor- 
age by  head  node)  and  assigns  as  much  flow  as  possible  to  its 
least  cost  leaving  (entering)  arc.   The  procedure  stops  when  an 
exhaustive  pass  of  the  nodes  makes  no  additional  flow  assignments. 
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The  starting  solution  achieved  is  not  necessarily  leasible. 
(F.g.,  "exhaustive  pass  sequential  source  minimum  start"  [Its].) 

Artificial  arcs  are  driven  from  the  basis  in  all  experiments 
using  a  big-M  method  (e.g.,  [b]).   This  choice  is  principally 
motivated  by  the  comparability  of  competitive  tests  between  pure 
and  generalized  network  codes  on  pure  and  generalized  network 
problems.    (A  two-phase  method  is  employed  in  production  use.) 

Choosing  the  best  big-M  value  is  a  bit  tricky.   The  smallest 
big-M  value  which  yields  a  feasible  optimal  solution  (if  one 
exists)  is  best  in  our  experience.   Small  big-M  values  may  tail 
to  produce  feasibility,  and  large  values  inflict  numerical 
difficulties.   In  practice,  a  default  value  is  used  and  an 
automatic  restart  recovery  is  applied  it  an  mteasible  solution 
persists.   If  a  restart  with  a  higher  big-M  value  tails  to 
reduce  the  total  inf eas ibility ,  a  terminally  inteasible  solution 
is  declared.   Figures  11  and  12  indicate  the  multiple  ot  maximum 
absolute  arc  cost  used  for  big-M  in  each  problem. 

Computational  tests  have  been  performed  on  various  computer 
systems.   The  times  reported  here  are  accurate  to  the  precision 
displayed  tor  IbM  37u/lbb-3  using  the  FORTRAN  H  compiler  with 
OPT(z).   Solution  times  exclude  input/output  overhead.   Ch'NNtT 
uses  double  precision  (IbM  RhJAL*b)  arithmetic  for  floating  point 
operations  and  storage. 

solution  times  are  given  in  Figure  11  for  the  pure  network 
test  problems.   Performance  is  given  tor  three  pure  network  codes 
(two  versions  ot  GNfcT  and  SUPFkk)  as  well  as  tor  several  repre- 
sentative generalized  network  prototype  systems.   we  can  thus 
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compare  the  best  generalized  network  scheme  with  the  best  pure 
network  code  exhibiting  equivalent  features  (GNfcT,  or  its 
aggregated  successor  variant  [6J). 

The  times  for  SUPEKK  also  provide  the  only  available  objec- 
tive means  lor  comparison  of  our  implementations  to  that  of 
Glover,  Hultz,  Klingman  and  btutz  [14];  they  report  that  their 
generalized  network  code,  N£TG,  is  about  as  fast  as  stiPfcKK  (a 
fast  out-ot-kilter  code  for  pure  networks  lij)  when  both  are  used 
to  solve  pure  networks.   Using  their  version  of  SUPtKK  on  our 
computer  we  have  shown  that  our  implementation  of  NETG  (called 
TLA  in  our  nomenclature)  is  at  least  as  efficient  as  their  claims 
for  NfcTG. 
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However,  in  these  tests  TLA  is  substantially  outperformed  by 
the  alternate  systems.   For  the  pure  network  problems,  best 
performance  is  achieved  by  (Figure  11): 

Contiguous  arcs  by  head  node, 

Candidate  queue  pricing, 

Preorder  traversal,  and 

Aggregated  successors. 
TLA  does  not  incorporate  any  of  these  features,  and  is  generally 
less  than  half  as  fast  as  competitors. 

Although  GENNET  (HgPX)  should  in  theory  rival  GNET  (HyPX) 
with  pure  network  problems,  the  overhead  of  testing  in  GENNET 
for  more  general  basis  structure  and  the  additional  computa- 
tional burden  of  floating  point  arithmetic  exact  a  performance 
penalty  of  about  15  percent. 

Figure  12  shows  solution  times  for  the  generalized  trans- 
shipment network  problems.   Note  that  the  starting  strategy  helps 
candidate  list  performance  and  hinders  the  candidate  queue. 
Arranging  arcs  contiguously  by  head  node  dominates  both  tail 
node  and  explicit  arc  list  designs.   The  candidate  queue  pro- 
vides good  performance  _if_  accompanied  with  contiguous  arcs  by 
head  node.   Preorder  traversal  continues  to  provide  better 
performance  than  triple  label  representation  in  all  design 
contexts.   Aggregated  successors  offer  a  pronounced  advantage. 

GENNET  (HyPX)  provides  best  overall  performance.   It 
offers  a  decided  advantage  on  problems  with  many  more  sinks 
than  sources,  a  situation  common  in  real  life. 
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Figure  13  displays  performance  of  GENNET(GN,  HPQX)  applied 
to  a  set  of  (GN)  problems  extracted  from  a  collection  of  real- 
world  LP/MIP  models  [1UJ.   Despite  the  slight  additional  floating 
point  arithmetic  required  to  solve  (GN),  GENNET  solves  these 
problems  much  faster  than  would  be  predicted  by  experience  with 
randomly  generated  GT  problems.   Tuning  of  the  pricing  mechanism 
greatly  enhances  this  difference. 


Problem 

Node-; 

Arcs 

Seconds 

Pivots 

AIRLP 

17U 

3,040 

2.62 

420 

CUAL 

170 

3,923 

1.80 

471 

STEEL 

422 

1,279 

.39 

499 

FOAM 

951 

4,953 

3.74 

1,258 

ODSAS(GN) 

1,431 

4,615 

3.22 

1,427 

ALUMINUM(GN) 

2,178 

7,216 

3.57 

2,794 

REFINE (GN) 

3,110 

6,617 

4.72 

3,322 

FOOD(GN) 

3,716 

13,907 

12.11 

7,004 

(Big-M  =  10  x  largest  cost  coefficient) 


Figure  13.   (GN)  Test  Problems 

(GN  rows  extracted  from  real-world  LP/MIP  models  [10 

with  null  columns  deleted  and  slack  arcs  added.) 
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Close  scrutiny  of  solution  trajectories  lends  some  insight 
into  GENNET's  good  performance.   GENNET  has  a  one-pass  itera- 
tion unless  a  cycle  is  formed;  a  cycle  is  formed  on  only  b-to-24 
percent  of  all  iterations  for  these  problem  sets — 5-to-lU  percent 
for  most  problems.   Also,  the  explicit  (non-aggregated)  subset  of 
the  nodes  is  remarkably  small,  seldom  numbering  much  more  than 
the  number  of  source  nodes.   Finally,  the  length  of  backpaths  is 
quite  short,  averaging  about  the  number  of  echelons  (path  length 
from  sources  to  sinks)  in  the  model,  or  just  more  than  2    in  these 
problems . 


56 


CONCLUSION 


The  generalized  network,  system  GENNh'T  is  small,  fast  and 
easy  to  modify.   Adaptations  have  already  included  using  GENNfcT 
in  a  system  to  solve  generalized  networks  with  complicating  side 
constraints  an/or  complicating  variables  (Mcbride  [3UJ).   GENNET 
has  also  been  incorporated  in  a  powerful  microcomputer-based 
network  optimization  system  by  Brown,  Duff  and  t'inley  [7,lbJ 
using  an  APPLE-II  host  and  PASCAL  implementation  language. 
Modifications  for  mixed  integer  generalized  networks  have  also 
been  tested  (though  not  with  care  sufficient  to  warrant 
publication  at  this  time).   GENNET  has  proven  to  be  a  worthy 
successor  of  GNET  [6]. 

Preorder  traversal  is  appealing  for  its  mathematical  and 
implementation  elegance,  and  has  proven  to  be  efficient  and 
flexible  for  generalized  networks  (as  it  was  tor  pure  networks). 
(Adolphson  and  Heum  [1]  have  also  suspected  this  and  have  in- 
dependently pursued  this  avenue.) 

Experience  shows  that  the  GENNET  design  pertorms  much  more 
efficiently  on  real  models  than  on  randomly  generated  test 
problems  of  nominally  equivalent  size;  this  design  is  also  tech- 
nically and  philosophically  compatible  with  the  various  systems 
we  have  devised  tor  solviny  other  more  general  classes  ot 
optimization  models. 
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The  EQRTRAN  programs  GENNET — (GT)  and  (GN)  versions  — 
(Copyright  1984)  are  licensed  to  researchers  lor  a  nominal 
charge  on  an  exclusive  use  basis.   For  further  information  write 
the  authors  via  P.O.  Box  1832,  Alexandria,  Virginia,  22313,  USA. 
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